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1. Introduction 


This article discusses some aspects of spectral methods for the solution 
of initial boundary value problems* The model problem can be formulated in 
the following way 


3U 

*5t 


- GU = 0 


U(x,0) = U®(x) 


( 1 . 1 ) 


where for each t, U(t) belongs to a Hilbert space H so that U(t) 

satisfies homogeneous boundary conditions and G is a linear spatial 
differential operator. There are three commonly used spectral methods for the 
space discretization of (1.1) - Galerkin, Tau and pseudospectral 

(collocation). Each one of three methods can be characterized by specifying a 
finite dimensional subspace Bjj c H and a projection operator Pj^ 

such that 

lim - uB = 0* (1*3) 

N- “ 

Using the operator Pjj results in a semidiscrete approximation of (1.1) 

^ U _ G U =0 
7t N N N 

Uj,(x,0) = uj(x) 


(1.4) 
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while 


“n " '■n " 


' "n ■ '■n “ 


S * 


(1.5) 


The commonly used basis functions of the subspace Bjj are related to 
Chebyshev or Legendre polynomials. 

Gjj is an operator from Bjj to Bj^; thus it can be viewed in the 
numerical procedure as an N x N matrix, the formal solution of (1.4) is 


Ujj(t) = exp(tGj^)U° . (1.6) 

When (1.4) is discretized in time by means of an explicit finite difference 
scheme, the time step is limited by a stability condition. It has been 

observed that the restriction on the time step At, for Chebyshev or Legendre 
method is of the form 

At = o(Iy-). 

8 

In fact when the equation (1.4) (for G “ discretized in space by the 

pseudospectral Chebyshev method and in time by the modified Euler scheme, then 
one encounters the stability condition [1] 


. (1.7) 

“ N 

The stability condition (1.7) is very stringent and has forced researchers to 
resort to implicit or semi-implicit time marching techniques thus complicating 
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the program and reducing the efficiency of the method. The stability 

condition (1.7) had been attributed to the well known CFL condition. Since 

the distribution of the grid points in any pseudospectral method is not 

uniform and “ o(n ^), then one should expect a stability condition of 

the form At ~ Ax . which agrees with (1.7). However, spectral methods are 
min 

global in nature since the solution at time step n+1 at a certain grid point 
depends on the solution at time step n at all the grid points. Therefore, 
an argument based on domain of dependence is not valid here. 

In this paper we analyze a pseudospectral method that does not have this 
severe limitation on the time step. This scheme is based on results obtained 
by M. Dubiner.^ In his paper, Dubiner has carried out a detailed analysis of 
the spectrum of the matrix % for the inflow problem 

(Vt ■ ^Vx " ° -1 < X < 1 (1.8) 

= U°(x) 

Uj^(l,t) = 0 

for several matrices Gjj resulting from various spectral approximations. He 
shows that if one uses the Tau method to solve (1.8) with Jacoby polynomials 
P' ’’^“'(x) as basis functions then the eigenvalues of G«, X., behave asymp- 
totically in the following way; 


^Dubiner, M., 1983, Tel-Aviv University, Tel-Aviv, Israel, personal 

communication. 
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= o(n^) o ^ 0 

x” = 0(N) a = 0 

1 


( 1 . 10 ) 


Using this result, we propose to show that it is possible to construct 
pseudospectral (collocation) algorithm for the solution of (1.8) such that the 
limitation on the time step is of the form 



( 1 . 11 ) 


It follows from Dubiner*s result (1.10) that in the Chebyshev case 

(a = ^ V2 ) > the domain in the complex plane which includes all the 

, 0 . 

eigenvalues of G^, has the size of order N . 
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While in the Legendre case (a = 0), the size of the domain Dj^ is of order 

N. 


Legendre Domain 



It is this difference in the size of the domains which results in different 
stability conditions. 

The choice of o = 0 in order to get (1.11), is because the inflow is 

from the right boundary.^ When the inflow is from the left boundary, one 
should choose an orthogonal polynomial with 3=0. For the case of inflow 
from both boundaries, the only choice is Legendre polynomial 

In Section 2 we derive a pseudospectral method that yields the same 
matrix G^j that corresponds to the Tau-Legendre method and proves the 
stability of the exact evolution operator exp(G^jt) . 

In Section 3 we analyze the solution of the fully discrete problem. Since 
spectral methods in space are highly accurate, it is desirable to have a 
similar accuracy in time as well. A scheme which has this property is 
explored in Section 4. And a slightly different approach for constructing 
Gji is described in Section 5. The algorithm based on this approach has some 
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advantages over the previous one from a programming point of view. On the 
other hand, instabilities occur when applied to a system of hyperbolic 
equations unless the boundary conditions are modified. 

In Section 6 we describe this phenomenon of instability and try to explain 
its origin. We also prove in this section that the first approach is stable 
without any modification of the boundary conditions. We conclude with Section 
7 giving some numerical results. 

2. The New Pseudospectral Method 

It has been shown [1] that when the Tau method is applied to the constant 
coefficient hyperbolic problem 


U(x,0) = U^(x) 

U(l,t) = 0 

then the numerical approximation satisfies exactly the equation 

■5^ - -5^ - x(t)Qjj(x) 

N 

” k=0 ^ ^ 


where 


Qj^(x) are any orthogonal polynomials. This has led many researchers to 
identify the Tau method with a collocation method where the collocation 
point Xj are the zeroes of Qj|(x) [2]. Note that the xj's lie in the 
open interval (-1,1). In order to construct this pseudospectral method, we 
define the following basis functions: 


Fjj(x) 

“ f;(x.Hx-x-T- 


1 < j < N 


( 2 . 1 ) 


where 


Fj^(x) = (x-Xj) •••(x-Xj^)(x-l) = (x-l)Qj^(x) 


( 2 . 2 ) 


It is easily verified that 


gj(Xk) = ®jk l< j, lc<N (2.3) 

and 

gj(l) = 0. 

Thus satisfies the right hand boundary condition. Using 

get the interpolation polynomial of the function U(x) 

N 

P„ U(x) = 5! U(x.)g.(x) (2;4) 

” j=l J J 


and its derivative 
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We now solve (1^8) by substituting Pjj U(x) instead of Ujj and satisfying 
the differential equation at the interior points xj. This, together with the 
homogeneous boundary condition (see remark at the end of the section in case 
of nonhomogeneous boundary conditions) results in the following set of 
equations: 

dU. N 

— = I g'.(xj^)U. l<k<N (2.6) 

j=l ■' 

where 

Uj = U(Xj). (2.7) 

(2.6) can be written in the matrix form 



( 2 . 8 ) 


(2.9) 



( 2 . 10 ) 


( 2 . 11 ) 
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The expression for (Gjj)ii is more involved. 


Define 


then for j = k we get 


R. (x) = x-x. 
1 1 


( 2 . 12 ) 


1 (x-x.)F'(x)-F (x) 

g-(xJ = fTCx T 2 


j Q^(x)+(x-l)Q^(x) 

F' (x.) x-x. 

N 1 x+x. j 

J 




J i=l 


N'“j' x>x.l* f 
Jf 1=1 


N N 




^ k=l i=l i=l 


N N 


^ k=l i=l 
k5^j i5tk 


N N 


B'V ii" |n 'i * EFl 'i >■ Fjr^ 1“ )f|niK«-i) E n 

^ J '^''j)i=l J k=l i=l ( ^ J k=l i=l 


N N 






k=l i=l 
k*3 i^^kjj 


N N 


k=l i=l 
k>^j i^kjj 


zn v«3> / 


(^N+l 


(2.13) 
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Thus, we finally have the following expressions for (Gjj 


(x,-UQJ(x.) X|^-x . 




'Vkj ■' 


1 


N+1 

lEx.-x. 

" J ^ 
i*3 


j=k 


(2.14) 


From the theory of the zeroes of Jacobi polynomials [6], we have the following 
identity 


E 

4=1 


1 


X. 

J 


-X. 

1 


+ 




0+1 

?Tx-=I) 


+ 


3+1 

2(xj+l) 


= 0 


2.15) 


where a and 3 are the powers in the expression of the weight function 
OL 6 

w(x) = (1-x) (1+x)^ of the Jacobi polynomials. In the Legendre case we have 
0=3=0; thus expression (2.14) can be simplified 




(yi)QN(Xk) 1 

(xj-l)Qj(x.) xj^-x . 


1 


x?-l 

J 


j^k 


j=k 


(2.16) 
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From the programming point of view it is convenient to define 


then 


where 


‘Vij ' 





j*k 


j=k 


(2.17) 


(2.18) 


(2.19) 


In order to use the operator Gjj, one has to store only two vectors - x,w 
where 

= \ 


(w) 


k 




( 2 . 20 ) 


The number of multiplications is 

N + + N = + 2N ~ . 


( 2 . 21 ) 


This number should be compared with CNlogN in the Chebyshev case (using 
FFT). For small N (up to N = 64) the two results are of the same order. 

The stability proof for the solution of the semidiscrete problem (1.8) is 
straightforward. Define the following vector norm 


■7 II 

w 



2 
V . 

J 


1 

T 


( 2 . 22 ) 
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where wj are the weights used in the Gause-Legendre quadrature, namely 


then 


2 [ q ^U )]2 

J=1 J=1 J J 


(2.23) 


(2.24) 


3 

The function U^(x)*g^ Uj^(x) is a polynomial of degree 2N-1 • Therefore, the 
Gauss-Legendre quadrature yields the exact value of the integral. Thus, we 
get 


7t 2 /!"» ["NJii = ”n“> - 

J=1 -1 


Using the boundary condition, " Oj results in 


Since 


we finally have 


k "^n"w < 0- 


= Bexp(G„t)u;jl 


N'"' N w 


lexp(Gj^t) 11^ < 1 


(2.26) 


(2.27) 


However, Dubiner's paper doesn't carry out a detailed analysis for the other 
two typical problems: 1) outflow, 2) inflow from both boundaries. It 
demonstrates how this analysis can be done and that the results concerning the 
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asymptotic behavior of the eigenvalues will be similar to (1.10). 
like to show how we define the operators Pjj and Gj^ for these two 
(1) Outflow 

A model problem for outlfow in both boundaries is 


(Vt ^ ° -1 < X < 1 

Uj^(x,0) = U°(x) 


This problem is well-posed without any boundary conditions. We 
define the basis functions as 


where 


gj(x) 


x-x. 


and 


consequently 


Qj^(x) = (x-Xj) •••(x-x^) (Legendre polynomial) 


N 

P U(x) = J] U(x .)g(x) . 

j=i J 




N 

I U(x.)g'.(x). 

j=l J J 


The elements of the matrix Gjj are 


We would 
problems. 

(2.28) 

therefore 

(2.29) 

(2.30) 

(2.31) 

(2.32) 


T 
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q;(x.) 




j^k 


x^-l 

J 


j=k 


Using the similarity transformation 


we get 


and 


<=N = "» S “n’ 


<Vlj - 


X .-X- 

J k 




2 
X . 

J 


x^-l 

J 


j*k 


j=k 


(2) Inflow From Both Boundaries 


The semidiscrete representation of the P.D.E. is 


(Vt- = ° 

Uj,(x,0) = U°(x) 

Ujj(-l,t) = Uj^(l,t) = 0 


(2.33) 


(2.34) 


(2.35) 


(2.36) 


(2.37) 
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Since the basis functions have to satisfy 


g.(l) = g.(-l) = 0 

we define 

, , 1 


where 


Sj^(x) = (x-Xq)(x-Xj) •••(x-Xj)(x-Xjj^j) = (x -l)Qjj(x) 


and 



N+1 


= 1 ) 


Pj, U(x) 


N 

I U(x.)g,(x). 

j=l J J 


The elements of the matrix Gj^ in this case are 


(Vjk 



j^k 


j=k 


6^ is similar to 


( 2 . 38 ) 


( 2 . 39 ) 


( 2 . 40 ) 


( 2 . 41 ) 


I 
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while 


Gm = G„ h"^ 
N N N N 




(2.43) 




^S>jk = 


(2.44) 


Remark ; When the boundary conditions are inhomogeneous, we have to modify 
our representation in the following way; For the right inflow problem we add 
another basis function 


_ 1 

“ F^l) 

N 


and we have 


(2.45) 


Sn+i^i) 


1 ; ^N+l^^k^ “ 1 £ k _< N 


(2.46) 


Thus, instead of (2.8) we get 


3t "» ■ “n "n * 


(2.47) 


f(t) = U(l,t) is the boundary conditon 


(2.48) 
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(W is defined by (2.20)) • When we have boundary conditions on both sides of 
the interval, we add two basis functions 




S„(x) j 

sJjTTT x+i 




Vx) j 

sjm x-1 


(2.50) 


hence 


gpC-l) =1 ; ^ ^ 1 ^ k _< N+1 (xjj^j=l) 

(2.51) 

gN+iU) = 1 ; gN+l^^k^ =0 0 < k < N (Xq= -1) 


Thus, instead of (2.8) we get 


(2.52) 

while 

f(t) = U(-l,t) is the left boundary condition 
g(t) = U(l,t) is the right boundary condition 

and 

^'^N^k “ (2.53) 




(2.54) 
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3. The Fully Discrete Solution 
The exact solution of (1.4) is 

\(x,t) = exp(tGj^)U° . (3.1) 

In [4] it has been shown that any explicit time algorithm can be represented 
as a polynomial approximation of the exact evolution operator exp(tGjj); thus 
the fully discrete solution of (1.1) is 

= ytG^)uO (3.2) 

where Hj|(z) is a polynomial of degree M which converges to e^ in the 
domain that includes all the eigenvalues of the matrix tGjj. The eigenvalues 
of tGjj are distinct^; therefore the corresponding eigenvectors are linearly 
independent and we can define a matrix Sjj whose columns are the eigenvectors 
of Gjj such that 

\ < 3 - 3 ) 

where Ejj is a diagonal matrix whose elements are 

® 

Therefore, if 

le^-V^)! 0 zei (3.5) 
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while Ijj is the domain in the complex plane which includes all the 
eignevalues of tGjj, then 


«U„ - V*!l > 0. (3.6) 

The relation between M and N depends mainly on three factors. 

1. The rate of convergence of Hj^(z) to e^. 

2. The size of the domain 1 ^. 

3. The norm of the matrices S.., . 

N N 

In [4] we find that for periodic problems where ^ ^ (the 

eigenvectors are orthogonal) one has to choose M such that the scalar 
function Hj^(z) resolves the exponent function e^ for In the case 

of boundary value problems, the analysis is much more complicated since the 
eignevectors are not orthogonal. We were not able to get an expression for 
the norms of 3 ^ and • However, numerical experiments verify the 

assumption that asymptotically one can get a sufficient condition relating 
M and N by carrying out an analysis based on the concept of resolution. 

Consider for example the modified Euler scheme. In the constant 
coefficients case, it is equivalent to the second order Taylor series method 

= (l+AtGj^+|(At)2 G^)v^ (3.7) 

or 

= (l+AtGj^+^(At)2 G^)^ U° . (3.8) 


If n is the number of time steps required to march to time level t, i.e., 
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At = t/n 


(3.9) 


V„(t) 


zn 


Thus, the numerical evolution operator is 


2n 


(M = 2n); (3.10) 


upon substituting z for Gj^t we get 


H^(z) = (l+ ^ z+ 

2n 


(3.11) 


Since 


- - [1* f * <3-12) 




0 < 0 < 1 


(3.13) 


e^-H^(z) I = n(Hj^(z)) R + low order terms (3.14) 


substituting (exp[-^)-R) for [h^j(z) ]” results in the following expression 
for the relative error E 


exp(z(n-l) /n) 


(for 1^1<1). (3.15) 


Using (3.13) in (3.15) gives 
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E 


T\exp((9-D^). 

n 


Thus, resolution of by Hjj(z) is achieved when 



(3.16) 


(3.17) 


(The magnitude of e is problem dependent). From (1.10), (3.9) and (3.10) it 
follows that in order to satisfy (3.17) we have to choose M such that 


or, equivalently 


M = 0(N^^^) 
At = 


(3.18) 


(3.19) 


The power 3/2 is due to the fact that the modified Euler scheme is second 
order accurate in time. A similar analysis for any explicit scheme which is 
p order accurate in time will yield the following condition 

P+1 

M = o(n P ). (3.20) 


It is obvious from this expression that using a scheme which has high accuracy 
in time will lead to the desired condition 

M = 0(N) ; (3.21) 


such a scheme is described in the next section 
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Remark ; Since we assume that resolution implies stability, conditions 
(3*18)-(3.10) are sufficient but not necessary# It is possible to get stable 
results while M satisfies M = 0(N) even for second order in time scheme as 
shown by the numerical results presented in Section 7. 


4* Highly Accurate Time Discretization 

The formal representation of an explicit fulldiscrete solution of (1#1) is 
given by (3#2). Since spectral methods in space are highly accurate, we would 
like to find a polynomial Hj^(z) that will yield high accuracy in time as 
well. Such a polynomial is described in [4] for pure initial value 

problems. It is based on approximating the function e^ by orthogonal 
polynomials. We would like to show how to implement this approach in the case 
of inflow - outflow boundary conditions. 

The main difference between the present case and the periodic one is the 
topological structure of the domain that includes the eigenvalues of 
In the periodic case we have (see remark at the end of the section) 

iReXjl < Cj ; ll^ X”| < C^CN) (4.1) 

N 

where X^^ are the eigenvalues of tGjj. 

(The constant Cj does not depend on N, while C 2 (N) = 0(N).) Vfhereas, in 
the boundary value case the eigenvalues satisfy^ 


iReX^l < Cj(N) ; \l^ xjl < C2(N) 


(4.2) 
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(Usually Cj^ 2^^^ “ O(N^). In Section 2 we have defined projection operators 
such that the related eigenvalues satisfy “ 0(N)). 

Accounting for this fact we have to modify the O.P.S. (Orthogonal 
Polynomials Scheme). 

Define 

S = max|Re(X^?)| (4.3) 

• X 

1 


R = maxi I (X^)|. 
. ‘ 1 ^* 

1 


(4.4) 


Since resolution of e^ by Hj^(z) means 


E 


max 

z 


e^-Hj^(z) 


z 

e 


< e 


z€D (domain of the eigenvalues) 


(4.5) 


for small enough e, we would like to choose Hj^(z) such that E is small 
for given M. Approximation based on the polynomials 4>|^(w) defined in [4] 
(i.e., orthogonal polynomials on the imaginary axis) will converge in D but 
will result in relatively large error E. Accounting for the fact that the 
denominator of E achieves its minimum in the left side of D, it is 
advisable to use a set of polynomials which are orthogonal on the line 
Re(z) = -s. Using this set of polynomials is equivalent to approximating 
e^ through the following change of variables 


Thus 


z = z + s. 


K(|) 


z -s z+s "S z -s 
e=e e =e e=e e 


(4.6) 


(4.7) 
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Therefore 


where 


= I b. (w) ; w = z/R 

k=0 

{ 1 k = 0 

2 k > 1 


(4.8) 

(4.9) 

(4.10) 


Jj^(R) is Bessel functions of order k. 4>^(w) satisfy the following 
recurrence relation 


<|>j^(w) = 2w(|>j^_j(w) + <l'j^_2(w) 

(4.11) 

<|>q(w) =1 ; <J>j(w) = w. 

Thus, substituting the operator 

( 4 . 12 ) 


instead of w in (4.8) results in the following approximation of the 
evolution operator 

- H„(W„) - ! (T^) (4.13) 

k=0 


and the fully discrete numerical solution is 
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= 

N 


M 

Ib 

|k=0 


k “r 


(4.14) 


where 


Numerical experiments show that while using the pseudospectral projection 
operator defined in Section 2 for the solution of the problems: 1) outflow, 
2) inflow from both boundaries, we have 


S » R. (4.15) 

Since, in this case the eigenvalues are grouped close to the real axis, the 
scheme described in [5] (for parabolic problems) will perform better than 
(4.14). (This conclusion is valid for any problem where we have an a priori 
information about the domain D similar to (4.15).) Hence we approximate the 
evolution operator exp(tGjj) in the following way: 

eKp(W„) . I^(G^) (4.16) 


where 




(4.17) 



-26- 




(4.18) 


is modified Bessel function of order k and Chebyshev 

polynomial. The numerical solution at time level t is 


N 


i: ■‘k "S 


k=0 


(4.19) 


Tk(Gjj)U^ is computed by using the recurrence relation 


Tj^(x) = 2xTj^_j(x) - \_2(x) k > 2 


(4.20) 


Thus 


Tq(x) =1 ; Tj(x) = X. 




** ‘'N \ * 


(4.21) 


Remark ; (4.1) has been proven in [5] for the periodic problem 

Uj. - a(x)U^ = 0 
U(x,0) = U^(x) 


where a(x) = sin2x. Similar technique can be applied to prove (4.1) when 
a(x) is any second degree trigonometric polynomial. Numerical experiments 



verify the assumption that (4.1) is valid for the general case, when a(x) 
can be represented as a finite degree trigonometric polynomial. 


5. Modification of the Pseudospectral Method 

The operator Pjj defined in Section 2 leads to some complexity when the 
boundary conditions are inhomogeneous, as shown by (2.45)-(2.54) . It is 
possible to overcome this difficulty by using a slightly different operator 

^N* 

Define the basis functions 

Fn(x) 

F^(Xj)(x-Xj) 


0 < j < N+1 (5.1) 


= -1; 


N+1 


= 1 ) 


where 


V"> 


(x+l)(x-Xj) •••(x-Xjj)(x-l) 


(x2-1)qN(x) 


(5.2) 


then 


N+1 

U„(x) = P„ U(x) = J] U(x.)g.(x) (5.3) 

” ” isO J J 


is polynomial of order N+1 interpolating U at the points xj, 
j = 0,1 , • •• , N,N+1 . Its derivative is given by 
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By using this projection operator we get the matrix Djj (the numerical 
derivative opertor) whose elements are: 




jk 




2x. 
i 

x ^-1 

J 


^*3 




(5.5) 


The matrix Dj^ can be written as 


where 


“« - "n “n “k 


(5.6) 


X .-X, 

J k 


k 5^ j 


%>jk = 


2x. 
i 

x^-l 

J 




\ 


QJU) 


0 < j = K < N+1 


j = k = 0 


j = k = N+1 


(5.7) 


and is a diagonal matrix 
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Thus 



k = 0 


0 < k < N+1 

2Q„(1) 

k = N+1 

■>» “n ' «N “n “n‘ 

®N 


(5.8) 


(5.9) 


where is a diagonal matrix 


^®N^kk 


1 

0 


0 _< k £ N 
k = N+1 


(5.10) 


Thus, we find that the algorithm is stable and the eignevalues of Gjg are 
<)• The main difference between the strategy used in Section 2 and the 
present one is the following: In the first case we follow exactly the P*D*E* 
and satisfy the equation only in the interior of the domain. The boundary 
conditions are satisfied by a proper choice of the basis functions. In the 
second case we satisfy the equation in the interior and boundary domain while 
imposing the boundary conditions at the end of each time step. Apparently, 
the first approach follows the P.D.E. more accurately than the second one. 
This statement will be made clear in the next section where we describe the 
solution of system of equations. 
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6 System of Hyperbolic Equations 
Consider the symmetric system 



U(-l) = f(t) ; U(l) = g(t) (6.1) 

0 .-C::) 

It is well known [3] that using a matrix of the type defined in the 

previous section in order to compute numerically the spatial derivatives leads 
to instability although the differential equations (5.1) are well-posed. It 
shows [3] how to stabilize the algorithm by adding numerical boundary 
conditions on the function V(x). This approach is based on the following 
argument. 

The characteristic variables are U+V and U-V and (6*1) is equivalent 
to the following diagonal set of equations ' 3 . 

(U-V)j. = - Y (U-V)^ (6.2) 

dx 3 • 

U+V is constant on the characteristic ~ ~ "2 constant on the 

characteristic ” T * Therefore, U+V should be given on the right 
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boundary and should be determined by the scheme on the left boundary. 
Similarly U-V should be given on the left boundary and should also be 

determined by the scheme on the right boundary. 

The scheme is stabilized by requiring that the values of U+V at 
X = -1 and U-V at x = 1 are not changed as a result of imposing the 
boundary conditions - U{1) and U(-l). 

It seems that this instability can be traced to the fact that by using the 
matrix we use the P.D.E. in the closed interval instead of the open 

interval as indicated by (6.1). By doing so we get errors on the boundaries 
for both U(x) and V(x). While the error in the function U(x) is 

immediately corrected by imposing the boundary conditions, the error in 
V(x) penetrates into the system through the characteristics and causes the 
instability. On the other hand, using the approach described in Section 2, we 
follow exactly the P.D.E. without imposing it on the boundaries; thus we do 
not expect this phenomenon of instability. This assumption is proved by the 
following theorem: 


Theorem: The solution of the semidiscrete problem 



-1 < X < 1 


(6.3) 
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discretized by the projection operator (2.38)-(2.40) for U and (2.29)-(2.31) 
for V, is stable. 


Proof Since (2 .38)-(2.40) Ujj is a polynomial of degree N+1, It can 
be represented as 

N+1^ 

Ujj = I (6.4) 

k=0 

where the P^^(x)'s are the normalized Legendre polynomials [Pj^(±l) = (±1)^], 
From (2.29)-(2.31) it follows that Vjj is a polynomial of degree N-1 and 
therefore 

N-1^ 

Vn = I vj^ Pk^^^ • (6.5) 

k=l 

Accounting for (6.4) and (6.5) the polynomials Ujj and Vj^ satisfy exactly 
the following equations 


<Vt 


¥V„ * <%>x * Vi<*> 


(6*6a) 






(6.6b) 


where Fjj are polynomials of degree N+1, N respectively, which vanish 

at the grid points. Therefore, we can write 


I would like to acknowledge ny advisor. Professor David Gottlieb, for this 
proof. 
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(ajX+bj^)Pj^(x) 
= ^2 


(6.7) 


( 6 . 8 ) 


In fact aj^, bj^ and a£ are given by 

„ - 2N+1 d " 

®i “n+r W “n+1 

bi = 4r '*N “ y(2N+l)uj^^.j 

A 

32 ~ “(2N+1 . 


(6.9a) 

(6.9b) 

(6.9c) 


(6.9) is proved by making use of the following relation satisfied by Legendre 
polynomials [6] 


xPj^(x) = 


(x) 

2N+1 N+r 


2N+1 




( 6 . 10 ) 


and the fact that ^ polynomial of degree N whose highest 

A 

coefficient is (2N+l)uj^^j . Thus, equating the coefficients of Pjj^j in 

(6.6a) results in (6.9a). Similarly, equating the coefficients of Pjj in 

(6.6a) results in (6.9b). Finally, equating the coefficients of Pj^ in 

(6.6b) results in (6.9c). 
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Define now the characteristic variables Rj^ and Sjj 


N+1 

k=0 


N+1^ 

Sn “ "n * - I fk<*> 

k=0 


N-1^ 

I T^. P^(x) + Uj, Pjj(x) 
k*0 


N-1^ 

I I*k(5c) + Ujj P^Cx) 

k=0 


(6*1 1 a) 


A 

*^N+l (6.11b) 


it is easily verified by (6.6), (6.7) that 

^^N^t "f^Vx (aiX+b^+a 2 )Pj^(x) (6.12a) 

^Vt = "T^Vx ^ (a^x+b^-a2)Pj^(x) (6.12b) 

multiplying (6.12b) by 3, adding it to (6.12a) and using the technique of 
equating highest coefficients of Pn- 1 results in 

^i ^ T5^'ar^®N-i'^^^N-i)* (6.13) 


Next, we use (6.12) to get 




1 

+ 3 / Rjj(ajX+bj-a 2 )Pjj dx. 


(6.14) 
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The first term on the R.H.S of (6.14) vanishes due to boundary conditions 
2 2 

(S^ = on the boundaries). 

From (6.9), (6*10), (6.11), (6.13) and the fact that Legendre polynomials 
are orthogonal, we get 


t4f “ 2 0^+1 4r 4+1 


(6.15) 


* 7 ^^N-1 ■ar^®N-i'^^’^N-i)^ * "ar 4 


where 



(6.16) 


On the other hand we have 


tIf / (4‘^^4)‘*^ "Tar 

-1 


TaJ:^3;2) 

k=0 


(6.17) 


4 Vl 


N+1 


Equating (6.15) and (6.17) gives us 


1 d 

T-ar 


N-2 


I a,(aj.3*r2) ♦ ^ 

k=0 


= 0 


(6.18) 


hence 


N-2 


k=0 


(6.19) 


where C is constant in time 
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From (6.11) and (6.19) we get 


N-2 

I a 

k=0 




.3"4) 


a 


N-l 


*2 

^N-1 


= C. 


( 6 . 20 ) 


All the terms on the L.H.S. of (6.20) are positive; therefore using (6.11) 
results in 


luj(t)| < Cl 
|vj(t)l < 


0 1 j 1 N-2 
0 <• j 1 N-l 


(6.21) 


while Cl, C2 are constant in time. Due to boundary conditions we have 


N+1 

U(-l) = I (-1)^ u^ = 0 
k=0 

N+1 

U(l) = I ° * 

k=0 

Thus for N even 


* “N - ■ - I 

k=0 


Un-1 + Uj, + Gjj^i = - I Uk 

k=0 


^N = - ^k 

k=0 

k even 


( 6 . 22 ) 


(6.23) 


(6.24) 
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N-3 

“N+1 * “n- 1 “ ~ “k • (6.25) 

k=l 
k odd 

• ^ A 

Since Uj^ for k = l,»»«,N-2 are bounded (6.21), then Ujj is bounded. In 

A 

order to show that (6.25) implies the boundedness of and we make 

use of (6.9a) and (6.12). Equating the two expressions for aj^ results in 

d " ^ N+1 d rr. _ 1 C ) rA oA^ 

ctF '^N+1 ~fl“ ^('‘n-1 Y ^N-iJ* (6.26) 

A 

^N-1 bounded; therefore in the limit we get 

4f “n+1 "" “n-1 N -*■ ») (6.27) 

A A 

(6*25) and (6.27) implies the boundedness of proof 

is concluded. Thus, we have proved stability of the semidiscrete problem 
(5.3), but (unfortunately) the domain of the eigenvalues of the related matrix 

A 

tGjq is proportional to and not to N as in the scalar case. This large 

domain will evidently result in a severe stability condition 

At = o(-ij-). 

N 


Hence, for the system case there is no difference between using Chebyshev or 
Legendre polynomials. 
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7* Numerical Results 

In this section we describe some numerical experiments whose results agree 
with the theory written in the previous sections. Throughout this section we 
use the following notations: 

N - Number of gridpoints (resolution in space). 

M - Degree of numerical evolution operator (resolution in time). 

X. - Eigenvalue of the related matrix tG^. 

The approximation in space is done by using the pseudospectral projection 
operator. The collocating points are the zeroes of the Nth degree 
polynomial. 

Table 1 presents the difference in the size of the domain of the 
eigenvalues while using Chebyshev or Legendre polynomials. The matrix tG^ 
whose eigenvalues we have computed is related to the problem 

- U =0 

t X 

U(x,0) = U°(x) (7.1) 

U(l,t) = f(t). 

We have taken t equal to 1 . 
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Table 1 


N 

Chebyshev 

max| 

Legendre 

max| 

8 

37.57 

7.0 

16 

150.0 

14.4 

32 

599.6 

29.8 


For the inflow problem 


- xU = 0 

t X 


U(x,0) = U®(x) 


(7.2) 


U(-l,t) = g(t); U(l,t) = f(t) 


the results were almost the same as in the previous table. For the third 
model problem of outflow 


+ xU =0 

t X 

U(x,0) = U°(x) 


(7.3) 


there is no difference between Chebyshev and Legendre* In both cases the 

eigenvalues are negative real numbers of order N, 

In Table 2 we compare the amount of work needed to achieve a certain 
degree of accuracy for Chebyshev and Legendre polynomials* The model problem 
is 


T 
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- U =0 

t X 


ttO/ ^ _ -X-1 

U (x) = e 


(7.4) 


U(l,t) = 0. 

The time marching technique is a fourth order Runge-Kutta. The solution is 
computed at time level T = 1.0. 


Table 2 


L 2 Error 

Chebyshev 

Legendre 


N 

M 

N 

M 

.32 X 10"2 

16 


240 

16 

48 

.246 X 10"3 

32 


960 

32 

100 


In Table 3 we carry out a similar comparison between Chebyshev and 
Legendre polynomials for the inflow problem 


- xU =0 

t X 

U(x,0) = exp(l/(x^-l)) (7.5) 


U(-l,t) = U(l,t) = 0 
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Table 3 


L 2 Error 

Chebyshev 

Legendre 


N 

M 

N 

M 

.645 X 10“^ 

16 


240 

16 


32 

.188 X 10"^ 

32 


960 

32 


64 


The next three tables are related to Section 4. The results presented here 
illustrate the high accuracy of the O.P.S. (Orthogonal Polynomials Scheme) 
compared to the second order Modified Euler scheme. Legendre polynomials are 
used for space approximations. In Table 4, the model problem is 


- U =0 

t X 


U°(x) = (x-l)^°, 


The solution is computed at t = 1.0. 


Table 4 


N 

Modified Euler 

O.P.S. 

M 

L 2 Error 

M 

L 2 Error 

16 

80 

0.1035 

36 

0.1660 X 10“5 

32 

160 

0.2388 X 10"^ 

72 

0.6836 X 10"^ 

64 

320 

0.5749 X 10"2 

144 

0.4247 X 10"^2 
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For the O.P.S., M was chosen such that the time error is of the same order as 
space error. This table shows clearly the overall spectral rate of 
convergence of the O.P.S. Comparing the O.P.S. to Modified Euler scheme with 
regard to the amount of work needed to achieve a certain degree of accuracy, 
one can use the fact that Modified Euler is second order in time. Thus, for 
N = 16 for example, an error of 0.1660 x 10”^ is achieved by the Modified 
Euler scheme when M satisfies 

M « 80(0.1035/0.1680 x lo"^)^ « 16000 
compared to 36 for O.P.S. 

The results in the next table are related to the inflow problem. 

- xU =0 

t X 

U(x,0) = (x^-1)^ 

U(-1,0) = U(1,0) = 0 
The solution is computed at t = 1.0. 
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Table 5 


N 

Modified Euler 

O.P.S 

M 

L2 Error 

M 

L2 Error 

16 

26 

0.1785 X 10"^ 

16 

0.1608 X 10"^ 

32 

52 

0.4265 X 10“2 

32 

0.1994 X 10”2 

64 

104 

0.8230 X 10“3 

64 

0.2004 X 10“^ 


Since the solution is very oscillatory, the advantage of using high order 
approximations (in space and time) is less significant than in the previous 
case. 

Table 6 describes the refinement procedure for the outflow problem. 


+ xU = 0 

t X 


U(x,0) = (x^-1)^ 


The solution is computed at t - 1.0. 


Table 6 


N 

Modified Euler 

O.P.S. 

M 

L2 Error 

M 

L2 Error 

16 

20 

0.9378 X 10“^ 

10 

0.8819 X 10"^ 

32 

40 

0.2463 X 10”3 

20 

0.2357 X 10"5 

64 

80 

0.6261 X 10"^ 

40 

0.2242 X 10“^° 
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8. Conclusion 

This paper has shown that it is possible to construct a pseudospectral 
method for initial boundary value scalar problems with stability condition 


it - o(^) 

rather than the familiar condition At = o(-^) . This improvement in the 

N 

stability condition does not hold in the case of system of equations or even 

in a scalar parabolic equation. Still, the fact of showing that using space 

discretization with Ax . = o(n 1 does not necessarily imply that 

min " 

At = o(n [for hyperbolic equation] gives us hope that there may be a way 
to overcome this drawback of using spectral methods for the numerical solution 
of nonperiodic problems. 
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